What is INLA? What is inlabru?

The short answer:

INLA is a fast method to do Bayesian inference with latent Gaussian models and inlabru is an R-package that implements this method with a flexible and simple interface.

The (much) longer answer:

cite some papers here

Where?

The software, information, examples and help can be found at

cite github for inlabru

• paper
• tutorials
• discussion group

So… Why should you use inlabru?

So… Why should you use inlabru?

• What type of problems can we solve?
• What type of models can we use?
• When can we use it?

To give proper answers to these questions, we need to start at the very beginning ..

The core

• We have observed something.

The core

• We have observed something.
• We have questions.

The core

• We have observed something.
• We have questions.
• We want answers!

How do we find answers?

We need to make choices:

• Bayesian or frequentist?
• How do we model the data?
• How do we compute the answer?

How do we find answers?

We need to make choices:

• Bayesian or frequentist?
• How do we model the data?
• How do we compute the answer?

These questions are not independent.

Bayesian or frequentist?

In this course we embrace the Bayesian perspective

• There are no “true but unknown” parameters !

Bayesian or frequentist?

In this course we embrace the Bayesian perspective

• There are no “true but unknown” parameters !
• Every parameter is described by a probability distribution!

Bayesian or frequentist?

In this course we embrace the Bayesian perspective

• There are no “true but unknown” parameters !
• Every parameter is described by a probability distribution!
• Evidence from the data is used to update the belief we had before observing the data!

Let’s formalize this a bit…

The elements of a inlabru friendly statistical model are:

1. The observational model \[ \begin{aligned} y_i|\eta_i, \sigma^2 & \sim\mathcal{N}(\eta_i,\sigma^2),\qquad i = 1,\dots,n\\ E(y_i|\eta_i, \sigma^2) & = \eta_i \end{aligned} \] Note: We assume that, given the linear predictor \(\eta\), the data are independent on each other! Data dependence is expressed through the components if the linear predictor.

Let’s formalize this a bit…

The elements of a inlabru friendly statistical model are:

1. The observational model \(y_i|\eta_i,\sigma^2\sim\mathcal{N}(\eta_i,\sigma^2),\qquad i = 1,\dots,n\)

2. A model for the linear predictor \[ E(y_i|\eta_i,\sigma^2) = \eta_i = \beta_0 + \beta_1x_i \]

Let’s formalize this a bit…

The elements of a inlabru friendly statistical model are:

1. The observational model \(y_i|\eta_i,\sigma^2\sim\mathcal{N}(\eta_i,\sigma^2),\qquad i = 1,\dots,n\)

2. A model for the linear predictor

\[ E(y_i|\eta_i,\sigma^2) = \eta_i = \color{red}{\boxed{\beta_0}} + \color{red}{\boxed{\beta_1x_i} } \]

Note 1: These are the components of our model! These explain the dependence structure of the data.

Let’s formalize this a bit…

The elements of a inlabru friendly statistical model are:

1. The observational model \(y_i|\eta_i,\sigma^2\sim\mathcal{N}(\eta_i,\sigma^2),\qquad i = 1,\dots,n\)

2. A model for the linear predictor \(\eta_i = \color{red}{\boxed{\beta_0}} + \color{red}{\boxed{\beta_1x_i} }\)

3. A prior for the model components \(\textbf{u}\) \[ \mathbf{u} = \{\beta_0, \beta_1\}\sim\mathcal{N}(0,\mathbf{Q}^{-1}) \] Note: These have always a Gaussian prior and are use to explain the dependencies among data!

Let’s formalize this a bit…

The elements of a inlabru friendly statistical model are:

1. The observational model \(y_i|\eta_i,\sigma^2\sim\mathcal{N}(\eta_i,\sigma^2),\qquad i = 1,\dots,n\)

2. A model for the linear predictor \(\eta_i = \color{red}{\boxed{\beta_0}} + \color{red}{\boxed{\beta_1x_i} }\)

3. A prior for the model components \(\mathbf{u} = \{\beta_0, \beta_1\}\sim\mathcal{N}(0,\mathbf{Q}^{-1})\)

4. A prior for the non-gaussian parameters \(\theta\) \[ \theta = \sigma^2 \]

Latent Gaussian Models (LGM)

Latent Gaussian Models (LGM)

Latent Gaussian Models (LGM)

Latent Gaussian Models (LGM)

\[ \pi(\mathbf{u},\theta|\mathbf{y})\propto \pi(\mathbf{y}|\mathbf{u},\theta)\pi(\mathbf{u}|\theta)\pi(\theta) \]

inlabru for linear regression

library(inlabru)
bodyfat[1:3,c("Abdomen","Bodyfat")]

  Abdomen Bodyfat
1    85.2    12.3
2    83.0     6.1
3    87.9    25.3

The code

# define model components
cmp =  ~ -1 + beta0(1) + beta1(Abdomen, model = "linear")

# define model predictor
eta = Bodyfat ~ beta0 + beta1

# build the observation model
lik = bru_obs(formula = eta,
              family = "gaussian",
              data = bodyfat)

# fit the model
fit = bru(cmp, lik)

inlabru for linear regression

bodyfat[1:3,c("Abdomen","Bodyfat")]

  Abdomen Bodyfat
1    85.2    12.3
2    83.0     6.1
3    87.9    25.3

The code

# define model components
cmp =  ~ -1 + beta0(1) + beta1(Abdomen, model = "linear")

# define model predictor
eta = Bodyfat ~ beta0 + beta1

# build the observation model
lik = bru_obs(formula = eta,
              family = "gaussian",
              data = bodyfat)

# fit the model
fit = bru(cmp, lik)

inlabru for linear regression

library(inlabru)
bodyfat[1:3,c("Abdomen","Bodyfat")]

  Abdomen Bodyfat
1    85.2    12.3
2    83.0     6.1
3    87.9    25.3

The code

# define model components
cmp =  ~ -1 + beta0(1) + beta1(Abdomen, model = "linear")

# define model predictor
eta = Bodyfat ~ beta0 + beta1

# build the observation model
lik = bru_obs(formula = eta,
              family = "gaussian",
              data = bodyfat)

# fit the model
fit = bru(cmp, lik)

inlabru for linear regression

library(inlabru)
bodyfat[1:3,c("Abdomen","Bodyfat")]

  Abdomen Bodyfat
1    85.2    12.3
2    83.0     6.1
3    87.9    25.3

The code

# define model components
cmp =  ~ -1 + beta0(1) + beta1(Abdomen, model = "linear")

# define model predictor
eta = Bodyfat ~ beta0 + beta1

# build the observation model
lik = bru_obs(formula = eta,
              family = "gaussian",
              data = bodyfat)

# fit the model
fit = bru(cmp, lik)

inlabru for linear regression

library(inlabru)
bodyfat[1:3,c("Abdomen","Bodyfat")]

  Abdomen Bodyfat
1    85.2    12.3
2    83.0     6.1
3    87.9    25.3

The code

# define model components
cmp =  ~ -1 + beta0(1) + beta1(Abdomen, model = "linear")

# define model predictor
eta = Bodyfat ~ beta0 + beta1

# build the observation model
lik = bru_obs(formula = eta,
              family = "gaussian",
              data = bodyfat)

# fit the model
fit = bru(cmp, lik)

inlabru for linear regression

Real-world datasets are more complicated!

Data can have several dependence structures: temporal, spatial,…

Using a Bayesian framework:

• Build (hierarchical) models to account for potentially complicated dependency structures in the data.

• Attribute uncertainty to model parameters and latent variables using priors.

Two main challenges:

• Need computationally efficient methods to calculate posteriors (this is where INLA helps!).
• Select priors in a sensible way (we’ll talk about this)

The second good news!

No matter how complicated is your model, the inlabru workflow is always the same 😃

# Define model components
comps <- component_1(...) + 
  component_2(...) + ...

# Define the model predictor
pred <- non_linear_function(component_1, 
                            component_2, ...)

# Build the observation model
lik <- bru_obs(formula = pred,
               family = ... ,
               data = ... ,
                ...)

# Fit the model
fit <- bru(comps, lik, ...)

The Tokyo rainfall data

One example with time series: Rainfall over 1 mm in the Tokyo area for each calendar day during two years (1983-84) are registered.

The Tokyo rainfall data

One example with time series: Rainfall over 1 mm in the Tokyo area for each calendar day during two years (1983-84) are registered.

The model

Stage 1 The observation model

\[ y_t|\eta_t\sim\text{Bin}(n_t, p_t),\qquad \eta_t = \text{logit}(p_t),\qquad i = 1,\dots,366 \] \[ n_{i} = \left\{ \begin{array}{lr} 1, & \text{for}\; 29\; \text{February}\\ 2, & \text{other days} \end{array}\right. \] \[ y_{i} = \begin{cases} \{0,1\}, & \text{for}\; 29\; \text{February}\\ \{0,1,2\}, & \text{other days} \end{cases} \]

• the likelihood has no hyperparameters

The model

Stage 1 The observation model

\[ y_t|\eta_t\sim\text{Bin}(n_t, p_t),\qquad \eta_t = \text{logit}(p_t),\qquad i = 1,\dots,366 \]

Stage 2 The latent field \[ \eta_t = \beta_0 + f(\text{time}_t) \]

• probability of rain on day \(i\) depends on \(x_i\)

• \(\beta_0\) is an intercept

• \(f(\text{time}_t)\) is a RW2 model (this is just a smoother). The smoothness is controlled by a hyperparameter \(\tau_f\)

The model

Stage 1 The observation model

\[ y_t|\eta_t\sim\text{Bin}(n_t, p_t),\qquad \eta_t = \text{logit}(p_t),\qquad i = 1,\dots,366 \]

Stage 2 The latent field \[ \eta_t = \beta_0 + f(\text{time}_t) \]

Stage 3 The hyperparameters

• The structured time effect is controlled by one parameter \(\tau_f\).

• We assign a prior to \(\tau_f\) to finalize the model.

inlabru for time series

Tokyo[1:3,]

  y n time
1 0 2    1
2 0 2    2
3 1 2    3

The code

# define model component
cmp =  ~ -1 + beta0(1) + time_effect(time, model = "rw2", cyclic = TRUE)

# define model predictor
eta = y ~ beta0 + time_effect

# build the observation model
lik = bru_obs(formula = eta,
              family = "binomial",
              Ntrials = n,
              data = Tokyo)

# fit the model
fit = bru(cmp, lik)

inlabru for time series

Tokyo[1:3,]

  y n time
1 0 2    1
2 0 2    2
3 1 2    3

The code

# define model component
cmp =  ~ -1 + beta0(1) + time_effect(time, model = "rw2", cyclic = TRUE)

# define model predictor
eta = y ~ beta0 + time_effect

# build the observation model
lik = bru_obs(formula = eta,
              family = "binomial",
              Ntrials = n,
              data = Tokyo)

# fit the model
fit = bru(cmp, lik)

inlabru for time series

Tokyo[1:3,]

  y n time
1 0 2    1
2 0 2    2
3 1 2    3

The code

# define model component
cmp =  ~ -1 + beta0(1) + time_effect(time, model = "rw2", cyclic = TRUE)

# define model predictor
eta = y ~ beta0 + time_effect

# build the observation model
lik = bru_obs(formula = eta,
              family = "binomial",
              Ntrials = n,
              data = Tokyo)

# fit the model
fit = bru(cmp, lik)

inlabru for time series

Tokyo[1:3,]

  y n time
1 0 2    1
2 0 2    2
3 1 2    3

The code

# define model component
cmp =  ~ -1 + beta0(1) + time_effect(time, model = "rw2", cyclic = TRUE)

# define model predictor
eta = y ~ beta0 + time_effect

# build the observation model
lik = bru_obs(formula = eta,
              family = "binomial",
              Ntrials = n,
              data = Tokyo)

# fit the model
fit = bru(cmp, lik)

inlabru for time series

Tokyo[1:3,]

  y n time
1 0 2    1
2 0 2    2
3 1 2    3

The code

# define model component
cmp =  ~ -1 + beta0(1) + time_effect(time, model = "rw2", cyclic = TRUE)

# define model predictor
eta = y ~ beta0 + time_effect

# build the observation model
lik = bru_obs(formula = eta,
              family = "binomial",
              Ntrials = n,
              data = Tokyo)

# fit the model
fit = bru(cmp, lik)

Example: disease mapping

We observed larynx cancer mortality counts for males in 544 district of Germany from 1986 to 1990 and want to make a model.

Take home message!